Subspace Identification of Modal Coordinate Time Series

The paper presents how direct estimation of the modal coordinate time series can be performed using a time domain subspace based system identification method. The method overcomes the traditional limitation on maximum number of modal coordinates to be estimated being less or equal to the applied number of sensors. This is achieved by utilizing the information on system order inherent in the state sequence applied as basis for the subspace identification of the state space system matrices. NOMENCLATURE K C M , , mass damping and stiffness matrices subscripts s or ,f refers to structure properties and force related properties respectively E D B A , , , state space system matrices, continuous time versions, subscripts a,v,q relates to acceleration, velocity and displacement ( ) ( ) ( ) t t t q q q & & & , , vector of generalised time dependent displacements, velocities and accelerations g K E D B A , , , , state space system matrices, discrete time, innovation form versions, including Kalman gain ( ) ( ) t t s , , , , , , , q q q f q q q f & & & & & & total dynamic forcing function, stochastic part and ( ) ( ) t t x x & , state vector, state vector time derivative ( ) ( ) t t r d f f , deterministic (measured) part, and residual (noise) part ( ) ( ) t t u y , system output and system input ( ) t d u B , input (load) influence matrix and deterministic inputs ( ) ( ) ( ) ( ) t t t k m w v e e , , , innovation and measurement noise, state output and process noise Y Y ˆ , matrix of output vectors, matrix of estimated k-step ahead predicted outputs X X ˆ , state sequence matrix, true and estimate U matrix of deterministic input vectors i W weighting matrix s L d L S S , coefficient matrices for deterministic and stochastic excitation in the extended state space model L O extended observability matrix 1 1 1 , , V S U matrices of an SVD G state estimation coefficient matrix R invertible scaling matrix t continuous time variable 0 t initial time k discrete time variable r m, number of measured outputs and inputs J L N , , number of samples, future horizon for identification, past horizon for instrumental variables 1 , − Ψ Ψ right and left eigenvector matrix j j ξ ψ , column vectors of right and left eigenvector matrices respectively Λ diagonal matrix of complex eigenvalues λ complex eigenvalue j ω is the damped circular natural frequency of mode j j 0 ω is the undamped circular natural frequency of mode j j α is the damping factor of mode j j ζ is the damping ratio of mode j ( ) 0 t j η is the initial complex modal coordinate of mode j corresponding to the eigenvalue pair λ λ j j , * and the initial state conditions x( ) t0 j φ is the initial modal phase of mode j i.e. ( ) ( ) 0 arg t a j j = φ T Super script T indicates transpose of a matrix jk θ is the phase of component k of right state space eigenvector j * Superscript * indicates complex conjugate ⊥ Super script ⊥ indicates orthogonal complement of a matrix { } o E expectation operator 1 − = i ô ^ over a parameter indicates estimate INTRODUCTION Modal analysis and system identification of vibrating structures have so far mainly dealt with the identification of natural frequencies, damping properties and mode shapes. However, when monitoring structures responding to natural excitation it is in many cases of great value to know the vibration amplitude for each excited mode. The standard straight forward way of achieving this is by decoupling the measured response using either experimental modeshapes or modeshapes obtained from a numerical model, as shown by e.g. Kaasen [1] or Hjelm et al [2]. However, one severe drawback with the methods presented so far is that only as many modal coordinates as there are sensors can be identified simultaneously. The present paper will introduce a method for determination of the modal coordinates where this restriction is lifted. The method which is based on the framework of subspace system identification makes it possible to simultaneously determine as many modal coordinates as there are identified natural frequencies and corresponding mode shapes in the data. The method was first applied for modal decomposition of the measured response of a drilling riser, Hoen and Moe [3]. However, in that paper details on the actual method where not given. In this paper we will give some details of the developed method for modal coordinate estimation. THE EQUATIONS OF MOTION FOR SYSTEM IDENTIFICATION OF VIBRATING STRUCTURES The dynamic response of a structural system can generally be modelled by a second order differential equation of dimension ( ) n n × as follows: ( ) ) , , , ( ) , , , ( ) ( ) ( ) ( t t t t t t s d s s s q q q f f q q q f q K q C q M & & & & & & & & & + = = + + (1) where q q & & & , and q are vectors of generalized acceleration, velocity and displacement, respectively. ) , , , ( t q q q f & & & is the forcing function which contains known (i.e. measured and thereby deterministic) excitation, ( ) t d f , and unknown (stochastic) excitation, ( ) t s , , , q q q f & & & . s s C M , and s K are the mass, damping and stiffness matrices of the structure. The stochastic part of the forcing function ) , , , ( t s q q q f & & & can be decomposed into a sum of elements being proportional to the acceleration, velocity and displacement respectively and a residual which contain all the other load components, also any non-linear effects: ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( t t t t t t t t r d s s s f f f f f q K q C q M q q q + + + + = + + & & & & & & (2) The first three elements of the right hand side of (2) are transferred to the left-hand side of the equation and expressed in terms of the acceleration, velocity and displacement respectively ( ) ( ) ( ) ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( t t t t t t t t r d f s f s f s f f q K K q C C q M M + = + + + + + & & & (3) ( ) t f M , ( ) t f C and ( ) t f K are mass (inertia), damping and stiffness effects caused by the external loading. In case of a structure submerged in water, the effects are known as hydrodynamic added mass and damping, and hydrostatic stiffness effects. The force related parts of the mass, damping and stiffness matrices are generally not time invariant. Therefore the general equation will contain time varying coefficient matrices. However, it is reasonable to assume that they may be regarded as approximately constant. This is at least reasonable in a time scale related to the time characteristics of the system, i.e. natural periods. Then we obtain the following secondorder differential equation ) ( ) ( ) ( ) ( ) ( t t t t t r d f u B Kq q C q M + = + + & & & (4) The mass matrix M is assumed positive definite. The damping matrix C may contain both viscous damping terms and gyroscopic terms. Gyroscopic terms may occur for e.g. risers with internal flow, and likewise for towed cables, se e.g. Blevins [4], and of course for rotating shafts etc. Thus, the damping matrix may be non-symmetric. The stiffness matrix K contains general stiffness properties. Normally the stiffness matrix will be symmetric. However, in certain flow-induced vibration problems, e.g. the classical flutter problem of airfoils, the equation of motion may be formulated to yield a non-symmetric stiffness matrix. d B is an input influence matrix characterising the locations and type of deterministic inputs ) (t u . The response of the dynamic system can be measured by e.g. accelerometers, inclinometers, rotation rate sensors, strain gages etc. A matrix output equation can thus be written as: ) ( ) ( ) ( ) ( ) ( t t t t t m q v a e q D q D q D y + + + = & & & (5) where the matrices v a D D , and q D are output influence matrices for acceleration, velocity and displacement respectively. ) (t m e is white measurement noise. The output influence matrices describe the relationship between the vectors q q q , , & & & and the measurement vector y . Thus, a measured output may be a combination of e.g. acceleration and rotation. This is in fact the case for accelerations measured with linear accelerometers mounted perpendicular to the longitudinal axis of a deep water riser as applied in offshore oil and gas exploration and production. For motions with a long period, the influence of the acceleration of gravity (the ” ( ) θ sin g⋅ ” component) may exceed the lateral acceleration in magnitude. This needs special attention during analysis of the measurements. In the case of interpreting system matrices identified or estimated from measured response, i.e. system identification, the system matrices cannot be assumed symmetric even if the tested system should yield symmetric matrices in theory. One major reason for non-symmetry in the identified matrices is that measurements always are imperfect and noisy. Thus, only under very special circumstances the eigenvalue problem of a system given by (4) will become symmetric and positive definite and thereby have real eigenvectors. In the general case complex eigenvectors occur. The eigenvalue problem corresponding to (4) can be solved in two ways, either by direct solution of the corresponding quadratic eigenvalue problem or as will be done here, by recasting (4) into a first order system in state space form. The state space model is a robust and good engineering model with a good numerical foundation for treating linear vibrating systems and it is as easy to understand as the second order approach. A STATE-SPACE MODEL Identification of the system parameters K C M , , , which in modal form are given by natural frequencies, modal damping ratios and mode shapes are not straightforward. The system identification methods applied in experimental modal analysis today are to a large extent based on a reformulation of the second order model (4) into a first order state-space description. See e.g. Juang [5]. In particular state space formulations have been applied for the purpose of system identification of offshore structures; see e.g. Hansteen [6], Hoen [7], Prevosto et al. [8]. Procedures for transformation of the second order model to statespace form can be found in textbooks on structural dynamics or system identification theory, see e.g. Hurty and Rubinstein [9], Juang [5] or Meirovitch [10]. By such procedures it is easy to see that it is always possible to represent a linear system given by (4) and (5) in state space form as follows: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) t t t t t t t t v Eu Dx y w Bu Ax x + + = + + = & (6) With reference to (4), (5) and (6) the following definitions apply ( ) ( ) ( )     = t t t q q x & is the state vector       − − = − − C M K M I 0 A 1 1 is the state transition matrix       − = − d B M 0 B 1 is the deterministic input matrix [ ] C M D D K M D D D 1 1 , − − − − = a v a q is the output matrix d a B M D E 1 − = is the deterministic feed-through matrix ( ) ( ) ( ) t t t m d a e w B M D v + = −1 is the state output noise ( ) ( )     − = − t t r f M 0 w 1 is the state process noise In case the state process noise ( ) t w is non-white, for practical purposes a state-space model can model the noise to yield a residual noise process that is white. This will add noise states to the state vector and corresponding terms to the matrices D B, A, . See e.g. Hoen [7] for details. Other forms of the state space representation are also possible depending on the definition of the state vector and the properties of the matrices C M, and K of (4). See e.g. Hurty and Rubinstein [9] or Laub and Arnold [11]. However, choice of formulation is only a matter of importance with respect to numerical implementation. They will all be related by simple coordinate transformations. A frequently applied alternative formulation to (6) in discrete time is the innovation form; see e.g. Ljung [12] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) k k k k k k k k g e u E x D y e K u B x A x + + = + + = +1 (7) where g K is the Kalman gain matrix and the innovation is defined as ( ) ( ) ( ) ( ) { } 1 E − − = k k k k y y y e where { } o E is the expectation operator. The system matrices E D B A , , , are the discrete time equivalents of the matrices E D B A , , , of (6). The innovation formulation is particularly useful for estimating the state vector time series, since it is known to yield optimal estimates of the state vector, see e.g. Maybeck [13]. THE STATE SPACE MODAL FORM The state space models (6) or (7) can be decoupled into a set of 2n uncoupled equations applying the eigenvalue decomposition of the state transition matrix, see Hoen [14, 3] for details. Λ AΨ Ψ = −1 (8) where Λ is the diagonal matrix of eigenvalues of A Ψ Ψ , 1 − is the left and right eigenvector matrices of A The eigenvector matrix of the state space model can be partitioned as       = Λ Ψ Ψ Ψ q q (9) where q Ψ is the components of the eigenvectors corresponding to the generalised displacements. Thus we obtain the following modal state space description by applying (8) to e.g. (6) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) t t t t t t t t v Eu η DΨ y w Ψ Bu Ψ Λη η + + = + + = − − 1 1 & (10) where ( ) ( ) t t x Ψ η 1 − = (11) is the complex vector of state modal coordinates. It is well known that the solutions to (6), (7) and (9) are composed of a homogeneous part associated with the initial conditions, and a steady state solution given by the future deterministic input and process noise. The solution to the homogeneous part is useful for interpretation of resonant vibrations such as e.g. lock-in Vortex Induced Vibrations of deep-water risers. The solution to the free vibration problem associated with (6) or (7) is known to be ( ) ( ) ( ) 0 0 t e t t t η Ψ x Λ − = (12) where ( ) 0 t η is a vector of complex coefficients or initial modal weights INTERPRETATION OF STATE SPACE MODAL RESPONSE The magnitude and the phase angle of the complex initial modal coordinate interpret as the initial modal amplitude and the initial modal phase angle. The elements of the initial modal coordinate vector ( ) 0 t η can therefore be written as ( ) ( ) 0 0 t t j j x ξ = η (13) where j ξ is the column vectors of the left eigenvector matrix 1 − Ψ . In matrix form the free vibration state response is given ( ) ( ) ( ) 0 0 1 , 0 t t t e t t t ≥ = − − Λ x Ψ Ψ x (14) Assume for simplicity of notation that Λ contains only complex eigenvalues, which then will appear in pairs as ( ) * , j j λ λ , where the asterisk denote complex conjugate. The free vibration response can then be expressed as the following sum over n components ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ∑ ∑